library(sp)
library(sf)
library(sfnetworks)
library(dodgr)
library(r5r)
library(tidygraph)
library(data.table)
library(dplyr)
library(purrr)
library(TSP)
library(openxlsx)
library(spatstat)
library(raster)
library(stars)
library(terra)
library(ggplot2)
library(ggspatial)
library(knitr)
library(tidyr)
library(stringr)
library(RColorBrewer)
library(foreach)
library(doFuture)
library(RColorBrewer)
ZAT <- st_read(dsn = "231031_ZAT_2023", layer = "231031_ZAT_2023") %>% st_transform(4326) %>% filter(UTAM != "N/A")
## Reading layer `231031_ZAT_2023' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota 2023\231031_ZAT_2023'
## using driver `ESRI Shapefile'
## Simple feature collection with 1215 features and 7 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY, XYZ
## Bounding box: xmin: 909732.5 ymin: 904238.5 xmax: 1113716 ymax: 1137247
## z_range: zmin: 0 zmax: 0
## Projected CRS: MAGNA-SIRGAS / Colombia Bogota zone
ZATCentroids <- st_centroid(ZAT)
setwd("C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota")
UTAM <- st_read(dsn = "Data", layer = "EMU2019_area")
## Reading layer `EMU2019_area' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 132 features and 57 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 569984 ymin: 492795.1 xmax: 626083.9 ymax: 557354.6
## Projected CRS: WGS 84 / UTM zone 18N
EF <- read.xlsx("Data/facteurs d'émissions.xlsx")
limites_bogota <- st_read(dsn = "Data", layer = "Localidad_Municipio_2017") %>% # Municipalities boundaries
st_transform(4326)
## Reading layer `Localidad_Municipio_2017' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 38 features and 6 fields
## Geometry type: MULTIPOLYGON
## Dimension: XY
## Bounding box: xmin: 564193.1 ymin: 483995.6 xmax: 634113.2 ymax: 570281.8
## Projected CRS: WGS 1984 UTM, Zone 18 North, Meter
setwd("C:/Users/hugot/Documents/Uniandes/EODH2023/EODH 2023 oficial/Publicacion1804/EODH/02_Base datos procesada/XLSX")
viajes <- read.xlsx("d. Modulo viajes.xlsx")
Hogares <- read.xlsx("a. Modulo hogares.xlsx")
Personas <- read.xlsx("c. Modulo personas.xlsx")
#Etapas <- read.xlsx("e. Modulo etapas.xlsx")
#Vehiculos <- read.xlsx("b. Modulo vehiculos.xlsx")
setwd("C:/Users/hugot/Documents/Uniandes/EODH2023")
#Esta base es necesaria para tener la ubicación del hogar
Hogares_loc <- read.xlsx(xlsxFile = "ANEXO 2_BASE DE DATOS EODH CON FACTORES DE EXPANSIÓN.xlsx", sheet = "hogar A-C-X-CF")
#Estas base son necesarias para deducir si el lugar de origen y destino de cada viaje es hogar u otro
Personas_viaj_no_procesada <- read.xlsx(xlsxFile = "ANEXO 2_BASE DE DATOS EODH CON FACTORES DE EXPANSIÓN.xlsx", sheet = "modulo D-E")
Viajes_no_procesada <- read.xlsx(xlsxFile = "ANEXO 2_BASE DE DATOS EODH CON FACTORES DE EXPANSIÓN.xlsx", sheet = "viajes")
setwd("C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota")
Routes <- st_read(dsn = "Data", layer = "Reseau_routier")
## Reading layer `Reseau_routier' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 100036 features and 12 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 565362.9 ymin: 486923.8 xmax: 633899.8 ymax: 569868
## Projected CRS: WGS 84 / UTM zone 18N
Routes <- st_cast(Routes, "LINESTRING") %>% st_transform(4326)
Transmi <- st_read(dsn = "Data", layer = "Transmilenio") %>% st_transform(4326)
## Reading layer `Transmilenio' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 13 features and 7 fields
## Geometry type: MULTILINESTRING
## Dimension: XY
## Bounding box: xmin: 588057.1 ymin: 503694.6 xmax: 606100.3 ymax: 527398.5
## Projected CRS: WGS 84 / UTM zone 18N
Transmi_stops <- st_read(dsn = "Data", layer = "estaciones-de-transmilenio") %>% st_transform(4326)
## Reading layer `estaciones-de-transmilenio' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota\Data'
## using driver `ESRI Shapefile'
## Simple feature collection with 140 features and 5 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: -74.20572 ymin: 4.55681 xmax: -74.04358 ymax: 4.768743
## Geodetic CRS: WGS 84
Se definen las categorías de modos que se usarán dentro del modelo de distancias
modos <- data.frame(modo_principal_desagrupado = unique(viajes$modo_principal_desagrupado),
modo_principal_typo = c("Auto", "Bicicleta", "A pie", "A pie", "Intermunicipal", "TransMilenio", "Transporte publico individual", "SITP Zonal", "Moto", "Alimentador", "Moto", "Auto", "Transporte publico individual", "Auto", "Transporte publico individual", "Bicitaxi", "SITP Zonal", "Transporte Escolar", "Auto", "Transporte informal", "Otro", "Otro", "Otro", "Transporte informal", "Auto", "Otro", "Bicicleta", "Cable", "Otro", "Patineta", "Bicicleta", "Otro", "Otro"),
modo_principal_comparado = c("Auto", "Bicicleta", "A pie", "A pie", "Transporte publico", "Transporte publico", "Taxi / Carro por aplicación", "Transporte publico", "Moto", "Transporte publico", "Moto", "Auto", "Taxi / Carro por aplicación", "Auto", "Taxi / Carro por aplicación", "Transporte informal", "Transporte publico", "Transporte Escolar", "Auto", "Transporte informal", "Otro", "Otro", "Otro", "Transporte informal", "Auto", "Otro", "Bicicleta", "Transporte publico", "Otro", "Otro", "Bicicleta", "Otro", "Otro"),
modo_principal_LATS = c("Private car", "Bike", "Walking", "Walking", "Regular bus", "BRT", "Taxi", "Regular bus", "Moto", "Regular bus", "Moto", "Private car", "Taxi", "Private car", "Taxi", "Paratransit", "BRT", "School bus", "Private car", "Paratransit", "Other", "Other", "Other", "Paratransit", "Private car", "Other", "Bike", "BRT", "Other", "Other", "Bike", "Other", "Other"))
viajes <- viajes %>% left_join(modos, by = "modo_principal_desagrupado")
Excluimos los viajes entrando o saliendo del area conformado por el Distrito de Bogotá y los 20 municipios que fueron encuestados. Para ello, se conservan solo los viajes que empiezan y terminan en las ZAT para los cuales tenemos un archivo SHP con su geometría y a los cuales les corresponde una UTAM. A la fecha, se usa el archivo con las 1,215 ZAT de 2023, y se reduce el número de viajes en la base de 100,174 a 96,054 (-4%).
#UTAM_presentes <- data.frame(UTAM = unique(c(Hogares$cod_utam_hg, viajes$utam_des, viajes$utam_ori))) %>% filter(UTAM != "No aplica")
#ZAT_presentes <- data.frame(ZAT = unique(c(Hogares$zat_hg, viajes$zat_des, viajes$zat_ori)))
#write.xlsx(UTAM_presentes, "UTAM_EODH_2023.xlsx")
#write.xlsx(ZAT_presentes, "ZAT_EODH_2023.xlsx")
viajes <- viajes[viajes$zat_des %in% ZAT$ZAT & viajes$zat_ori %in% ZAT$ZAT,]
Necesitamos recuperar las coordenadas de los hogares desde la versión de trabajo de la encuesta.
Hogares <- Hogares %>%
left_join(Hogares_loc[,c("i7a-Latitude", "i7a-Longitude", "KEY")], by = c("key_hg" = "KEY")) %>%
st_as_sf(coords = c("i7a-Longitude", "i7a-Latitude"), crs = 4326)
Hogares_loc <- Hogares_loc[,c("i7a-Latitude", "i7a-Longitude", "cod_dane","KEY")] %>%
st_as_sf(coords = c("i7a-Longitude", "i7a-Latitude"), crs = 4326)
#st_write(Hogares_loc, "Hogares_loc_mza.shp")
#Definición de la ventana
zoom_level <- 9.5
lon_span <- 360 / 2^zoom_level
lat_span <- 360 / 2^zoom_level
zoom_to_Bogota <- c(-74.12, 4.72) # Bogota
lon_bounds_Bogota <- c(zoom_to_Bogota[1] - lon_span / 2, zoom_to_Bogota[1] + lon_span / 2)
lat_bounds_Bogota <- c(zoom_to_Bogota[2] - lat_span / 2, zoom_to_Bogota[2] + lat_span / 2)
UTAM <- UTAM %>% st_transform(4326)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, fill = "white")+
geom_sf(data = Hogares_loc, col = "red", size = 0.5)+
labs(title = "Ubicación de los hogares de la muestra",
subtitle = "Bogotá DC + 20 municipios",
caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023",
col = "UPZ") +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(x = "", y = "") +
theme(legend.position = "right",
legend.text = element_text(size=10),
legend.title = element_text(size=15),
plot.title = element_text(size=15),
plot.caption = element_text(size = 7, face = "italic", hjust = 0, vjust = 12),
legend.text.align = 1)+
annotation_scale(location = "br", height = unit(0.12, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(0.8, "cm"))
Finalmente necesitamos identificar para cada viaje si sus puntos de origen y destino son el hogar de la persona u otro. Para hacerlo, necesitamos buscar en la versión de trabajo de la base de viajes la variable que contiene esta información para el lugar de destino únicamente (lastimosamente no la hay para el lugar de origen). Luego, para el lugar de origen del primer viaje del día, buscamos esta información en la version de trabajo de la base de viajes. Este procedimiento busca aumentar la precisión de los viajes que comienzan o terminan en el hogar
#Lugar destino
viajes <- viajes %>% left_join(Viajes_no_procesada[,c("d24","KEY")], by = c("key_viaje" = "KEY"))
colnames(viajes)[colnames(viajes) == "d24"] <- "lugar_destino"
viajes$lugar_destino[viajes$lugar_destino == 1] <- "Hogar"
viajes$lugar_destino[viajes$lugar_destino == 2] <- "Otro lugar"
#Lugar origen (funciona solo para el primer viaje del día)
Personas_viaj_no_procesada$KEY_2 <- gsub(pattern = "person_d_e", replacement = "B/person_b", Personas_viaj_no_procesada$KEY)
Personas <- Personas %>% left_join(Personas_viaj_no_procesada[,c("d7", "KEY_2")], by = c("key_persona" = "KEY_2"))
colnames(Personas)[colnames(Personas) == "d7"] <- "lugar_ini_dia"
Personas$lugar_ini_dia[Personas$lugar_ini_dia == 1] <- "Hogar"
viajes <- viajes %>% left_join(Personas[,c("cod_per", "lugar_ini_dia")], by = c("cod_pers" = "cod_per"))
viajes$lugar_origen <- "Otro lugar"
viajes$lugar_origen[viajes$lugar_ini_dia == "Hogar" & viajes$orden_vj == 1] <- "Hogar"
#Control (no debe haber viajes Hogar --> Hogar)
viajes$test_origen_destino[viajes$lugar_origen == "Hogar" & viajes$lugar_destino == "Hogar"] <- "HOGAR_HOGAR"
viajes$test_origen_destino[viajes$lugar_origen == "Hogar" & viajes$lugar_destino == "Otro lugar"] <- "HOGAR_Otro"
viajes$test_origen_destino[viajes$lugar_origen == "Otro lugar" & viajes$lugar_destino == "Hogar"] <- "Otro_HOGAR"
viajes$test_origen_destino[viajes$lugar_origen == "Otro lugar" & viajes$lugar_destino == "Otro lugar"] <- "Otro_Otro"
print(table(viajes$test_origen_destino))
##
## HOGAR_HOGAR HOGAR_Otro Otro_HOGAR Otro_Otro
## 1 45499 45721 4833
#Encontramos un (1) viaje "Hogar --> Hogar". Quitamos todos los viajes de la persona que lo realiza del registro de viajes
excl <- viajes$cod_pers[viajes$test_origen_destino == "HOGAR_HOGAR"]
viajes <- viajes[!(viajes$cod_pers %in% excl),]
# Ponderación de los segmentos y restricción al principal componente conexo
net_mot <- weight_streetnet(Routes, wt_profile = "motorcar", type_col = "type")
net_mot <- net_mot[net_mot$component == 1,]
net_cycl <- weight_streetnet(Routes, wt_profile = "bicycle", type_col = "type")
net_cycl <- net_cycl[net_cycl$component == 1,]
net_walk <- weight_streetnet(Routes, wt_profile = "foot", type_col = "type")
net_walk <- net_walk[net_walk$component == 1,]
Transmi <- st_cast(Transmi, "LINESTRING")
Transmi_stops <- Transmi_stops %>% st_transform(4326)
# Inicialización
net_transmi = as_sfnetwork(Transmi, directed = FALSE) %>%
st_transform(4326) %>%
activate("edges") %>%
mutate(weight = edge_length())
# Subdivisión de ejes
net_transmi = convert(net_transmi, to_spatial_subdivision)
# Inclusión (blend) de las estaciones a la red
net_transmi = st_network_blend(net_transmi, Transmi_stops)
# Cálculo de la nueva longitud de cada eje
net_transmi <- net_transmi %>% activate("edges") %>% mutate(weight = edge_length())
#Allocating 8 GB of memory to Java (2 GB is a minimum, adjust according to the hardware)
options(java.parameters = "-Xmx12G")
#gc()
#Initiating Java
#rJava::.jinit()
#rJava::.jcall("java.lang.System", "S", "getProperty", "java.version")
#Setting up the multimodal network. Requires a .zip archive containing the GTFS and a .pbf file with the road network.
#It is quite long the first time (~15 min), then it reloads the previously created network.
r5r_core <- setup_r5(data_path = "C:/Users/hugot/Documents/Rennes 2/Calcul distance SIG/Reprise perso/RMD Bogota/Data")
# set inputs
mode <- c("WALK", "TRANSIT")
departure_datetime <- as.POSIXct("13-05-2019 14:00:00",
format = "%d-%m-%Y %H:%M:%S")
#Putting coordinates in the appropriate format for r5r
HogaresR5R <- cbind(as.data.frame(st_coordinates(Hogares)), Hogares$cod_hog)
colnames(HogaresR5R) <- c("lon", "lat", "id")
ZATR5R <- cbind(as.data.frame(st_coordinates(ZATCentroids)), ZATCentroids$ZAT)
colnames(ZATR5R) <- c("lon", "lat", "id")
Se calcula la longitud de cada viaje de acuerdo a los siguiente:
| Modo principal | Red empleada |
|---|---|
| Transmilenio | Transmilenio |
| SITP Zonal | SITP |
| Alimentador | SITP |
| A pie | Red vial (ponderación peatón) |
| Bicicleta | Red vial (ponderación ciclista) |
| Cualquier otro | Red vial (ponderación automovilista) |
print(Sys.time())
#5 min
#LINEA RECTA (para sustituir NA)
# 1 --> 2 (Hogar --> ZAT)
mat_straight_12 <- as.data.frame(st_distance(x = st_geometry(Hogares), y = st_geometry(ZATCentroids)))
rownames(mat_straight_12) <- Hogares$cod_hog
colnames(mat_straight_12) <- ZAT$ZAT
mat_straight_12_table <- as.data.frame(as.table(as.matrix(mat_straight_12)))
colnames(mat_straight_12_table) <- c("cod_hog", "zat_des", "dist_straight")
# 2 --> 1 (ZAT --> Hogar)
mat_straight_21 <- as.data.frame(st_distance(x = st_geometry(ZATCentroids), y = st_geometry(Hogares)))
rownames(mat_straight_21) <- ZAT$ZAT
colnames(mat_straight_21) <- Hogares$cod_hog
mat_straight_21_table <- as.data.frame(as.table(as.matrix(mat_straight_21)))
colnames(mat_straight_21_table) <- c("zat_ori", "cod_hog", "dist_straight")
# 2 --> 2 (ZAT --> ZAT)
mat_straight_22 <- as.data.frame(st_distance(x = st_geometry(ZATCentroids), y = st_geometry(ZATCentroids)))
rownames(mat_straight_22) <- ZAT$ZAT
colnames(mat_straight_22) <- ZAT$ZAT
mat_straight_22_table <- as.data.frame(as.table(as.matrix(mat_straight_22)))
colnames(mat_straight_22_table) <- c("zat_ori", "zat_des", "dist_straight")
#TRANSMILENIO
# 1 --> 2 (Hogar --> ZAT)
mat_transmi_12 <- st_network_cost(net_transmi, from = st_geometry(Hogares), to = st_geometry(ZATCentroids))
rownames(mat_transmi_12) <- Hogares$cod_hog
colnames(mat_transmi_12) <- ZAT$ZAT
mat_transmi_12_table <- as.data.frame(as.table(mat_transmi_12))
colnames(mat_transmi_12_table) <- c("cod_hog", "zat_des", "dist_transmi")
# 2 --> 1 (ZAT --> Hogar)
mat_transmi_21 <- st_network_cost(net_transmi, from = st_geometry(ZATCentroids), to = st_geometry(Hogares))
rownames(mat_transmi_21) <- ZAT$ZAT
colnames(mat_transmi_21) <- Hogares$cod_hog
mat_transmi_21_table <- as.data.frame(as.table(mat_transmi_21))
colnames(mat_transmi_21_table) <- c("zat_ori", "cod_hog", "dist_transmi")
# 2 --> 2 (ZAT --> ZAT)
mat_transmi_22 <- st_network_cost(net_transmi, from = st_geometry(ZATCentroids), to = st_geometry(ZATCentroids))
rownames(mat_transmi_22) <- ZAT$ZAT
colnames(mat_transmi_22) <- ZAT$ZAT
mat_transmi_22_table <- as.data.frame(as.table(mat_transmi_22))
colnames(mat_transmi_22_table) <- c("zat_ori", "zat_des", "dist_transmi")
#MOTORIZADO
# 1 --> 2 (Hogar --> ZAT) (10 secs)
mat_mot_12 <- dodgr_dists(graph = net_mot, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_mot_12) <- Hogares$cod_hog
colnames(mat_mot_12) <- ZAT$ZAT
mat_mot_12[is.na(mat_mot_12)] <- mat_straight_12[is.na(mat_mot_12)]
mat_mot_12_table <- as.data.frame(as.table(mat_mot_12))
colnames(mat_mot_12_table) <- c("cod_hog", "zat_des", "dist_car")
# 2 --> 1 (ZAT --> Hogar) (10 secs)
mat_mot_21 <- dodgr_dists(graph = net_mot, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(Hogares)), shortest = TRUE)
rownames(mat_mot_21) <- ZAT$ZAT
colnames(mat_mot_21) <- Hogares$cod_hog
mat_mot_21[is.na(mat_mot_21)] <- mat_straight_21[is.na(mat_mot_21)]
mat_mot_21_table <- as.data.frame(as.table(mat_mot_21))
colnames(mat_mot_21_table) <- c("zat_ori", "cod_hog", "dist_car")
# 2 --> 2 (ZAT --> ZAT) (10 secs)
mat_mot_22 <- dodgr_dists(graph = net_mot, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_mot_22) <- ZAT$ZAT
colnames(mat_mot_22) <- ZAT$ZAT
mat_mot_22[is.na(mat_mot_22)] <- mat_straight_22[is.na(mat_mot_22)]
mat_mot_22_table <- as.data.frame(as.table(mat_mot_22))
colnames(mat_mot_22_table) <- c("zat_ori", "zat_des", "dist_car")
#BICI
# 1 --> 2 (Hogar --> ZAT) (10 secs)
mat_cycl_12 <- dodgr_dists(graph = net_cycl, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_cycl_12) <- Hogares$cod_hog
colnames(mat_cycl_12) <- ZAT$ZAT
mat_cycl_12[is.na(mat_cycl_12)] <- mat_straight_12[is.na(mat_cycl_12)]
mat_cycl_12_table <- as.data.frame(as.table(mat_cycl_12))
colnames(mat_cycl_12_table) <- c("cod_hog", "zat_des", "dist_bike")
# 2 --> 1 (ZAT --> Hogar) (10 secs)
mat_cycl_21 <- dodgr_dists(graph = net_cycl, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(Hogares)), shortest = TRUE)
rownames(mat_cycl_21) <- ZAT$ZAT
colnames(mat_cycl_21) <- Hogares$cod_hog
mat_cycl_21[is.na(mat_cycl_21)] <- mat_straight_21[is.na(mat_cycl_21)]
mat_cycl_21_table <- as.data.frame(as.table(mat_cycl_21))
colnames(mat_cycl_21_table) <- c("zat_ori", "cod_hog", "dist_bike")
# 2 --> 2 (ZAT --> ZAT) (10 secs)
mat_cycl_22 <- dodgr_dists(graph = net_cycl, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_cycl_22) <- ZAT$ZAT
colnames(mat_cycl_22) <- ZAT$ZAT
mat_cycl_22[is.na(mat_cycl_22)] <- mat_straight_22[is.na(mat_cycl_22)]
mat_cycl_22_table <- as.data.frame(as.table(mat_cycl_22))
colnames(mat_cycl_22_table) <- c("zat_ori", "zat_des", "dist_bike")
#WALK
# 1 --> 2 (Hogar --> ZAT) (10 secs)
mat_walk_12 <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_walk_12) <- Hogares$cod_hog
colnames(mat_walk_12) <- ZAT$ZAT
mat_walk_12[is.na(mat_walk_12)] <- mat_straight_12[is.na(mat_walk_12)]
mat_walk_12_table <- as.data.frame(as.table(mat_walk_12))
colnames(mat_walk_12_table) <- c("cod_hog", "zat_des", "dist_walk")
# 2 --> 1 (ZAT --> Hogar) (10 secs)
mat_walk_21 <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(Hogares)), shortest = TRUE)
rownames(mat_walk_21) <- ZAT$ZAT
colnames(mat_walk_21) <- Hogares$cod_hog
mat_walk_21[is.na(mat_walk_21)] <- mat_straight_21[is.na(mat_walk_21)]
mat_walk_21_table <- as.data.frame(as.table(mat_walk_21))
colnames(mat_walk_21_table) <- c("zat_ori", "cod_hog", "dist_walk")
# 2 --> 2 (ZAT --> ZAT) (10 secs)
mat_walk_22 <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(ZATCentroids)), shortest = TRUE)
rownames(mat_walk_22) <- ZAT$ZAT
colnames(mat_walk_22) <- ZAT$ZAT
mat_walk_22[is.na(mat_walk_22)] <- mat_straight_22[is.na(mat_walk_22)]
mat_walk_22_table <- as.data.frame(as.table(mat_walk_22))
colnames(mat_walk_22_table) <- c("zat_ori", "zat_des", "dist_walk")
print(Sys.time())
colnames(viajes)[colnames(viajes) == "cod_hg"] <- "cod_hog"
#Separando la base en 3 sub-bases de acuerdo al lugar de origen y destino
viajes_12 <- viajes %>% filter(test_origen_destino == "HOGAR_Otro")
viajes_21 <- viajes %>% filter(test_origen_destino == "Otro_HOGAR")
viajes_22 <- viajes %>% filter(test_origen_destino == "Otro_Otro")
#Unión por columnas múltiples (15 min)
print(Sys.time())
viajes_12 <- merge(viajes_12, mat_mot_12_table, by = c("cod_hog", "zat_des"))
viajes_21 <- merge(viajes_21, mat_mot_21_table, by = c("cod_hog", "zat_ori"))
viajes_22 <- merge(viajes_22, mat_mot_22_table, by = c("zat_ori", "zat_des"))
viajes_12 <- merge(viajes_12, mat_cycl_12_table, by = c("cod_hog", "zat_des"))
viajes_21 <- merge(viajes_21, mat_cycl_21_table, by = c("cod_hog", "zat_ori"))
viajes_22 <- merge(viajes_22, mat_cycl_22_table, by = c("zat_ori", "zat_des"))
viajes_12 <- merge(viajes_12, mat_walk_12_table, by = c("cod_hog", "zat_des"))
viajes_21 <- merge(viajes_21, mat_walk_21_table, by = c("cod_hog", "zat_ori"))
viajes_22 <- merge(viajes_22, mat_walk_22_table, by = c("zat_ori", "zat_des"))
viajes_12 <- merge(viajes_12, mat_transmi_12_table, by = c("cod_hog", "zat_des"))
viajes_21 <- merge(viajes_21, mat_transmi_21_table, by = c("cod_hog", "zat_ori"))
viajes_22 <- merge(viajes_22, mat_transmi_22_table, by = c("zat_ori", "zat_des"))
viajes_12 <- merge(viajes_12, mat_straight_12_table, by = c("cod_hog", "zat_des"))
viajes_21 <- merge(viajes_21, mat_straight_21_table, by = c("cod_hog", "zat_ori"))
viajes_22 <- merge(viajes_22, mat_straight_22_table, by = c("zat_ori", "zat_des"))
print(Sys.time())
## 27 HORAS
print(Sys.time())
#1-> 2 (Hogar -> ZAT)
viajes_12_sitp <- viajes_12[viajes_12$modo_principal_typo == "SITP Zonal",]
viajes_12_origen <- cbind(id = viajes_12_sitp$cod_hog, HogaresR5R[match(viajes_12_sitp$cod_hog, HogaresR5R$id), 1:2, drop = FALSE])
viajes_12_destino <- cbind(id = viajes_12_sitp$zat_des, ZATR5R[match(viajes_12_sitp$zat_des, ZATR5R$id), 1:2, drop = TRUE])
viajes_12_sitp$dist_sitp = 0
for(i in 1:nrow(viajes_12_sitp)){
#print(i)
det <- detailed_itineraries(r5r_core = r5r_core,
origins = viajes_12_origen[i,],
destinations = viajes_12_destino[i,],
mode = mode,
departure_datetime = departure_datetime,
max_walk_time = 180,
max_trip_duration = 300,
shortest_path = TRUE)
if(!is.na(det$total_distance[1])){
det <- det %>% filter(mode == "BUS")
viajes_12_sitp$dist_sitp[i] <- sum(det$distance)
}
}
print(Sys.time())
write.xlsx(viajes_12_sitp, "viajes_12_sitp.xlsx")
#2-> 1 (ZAT -> Hogar)
viajes_21_sitp <- viajes_21[viajes_21$modo_principal_typo == "SITP Zonal",]
viajes_21_origen <- cbind(id = viajes_21_sitp$zat_ori, ZATR5R[match(viajes_21_sitp$zat_ori, ZATR5R$id), 1:2, drop = TRUE])
viajes_21_destino <- cbind(id = viajes_21_sitp$cod_hog, HogaresR5R[match(viajes_21_sitp$cod_hog, HogaresR5R$id), 1:2, drop = FALSE])
viajes_21_sitp$dist_sitp = 0
for(i in 1:nrow(viajes_21_sitp)){
#print(i)
det <- detailed_itineraries(r5r_core = r5r_core,
origins = viajes_21_origen[i,],
destinations = viajes_21_destino[i,],
mode = mode,
departure_datetime = departure_datetime,
max_walk_time = 180,
max_trip_duration = 300,
shortest_path = TRUE)
if(!is.na(det$total_distance[1])){
det <- det %>% filter(mode == "BUS")
viajes_21_sitp$dist_sitp[i] <- sum(det$distance)
}
}
print(Sys.time())
write.xlsx(viajes_21_sitp, "viajes_21_sitp.xlsx")
#2-> 2 (ZAT -> ZAT)
viajes_22_sitp <- viajes_22[viajes_22$modo_principal_typo == "SITP Zonal",]
viajes_22_origen <- cbind(id = viajes_22_sitp$zat_ori, ZATR5R[match(viajes_22_sitp$zat_ori, ZATR5R$id), 1:2, drop = TRUE])
viajes_22_destino <- cbind(id = viajes_22_sitp$cod_hog, HogaresR5R[match(viajes_22_sitp$cod_hog, HogaresR5R$id), 1:2, drop = FALSE])
viajes_22_sitp$dist_sitp = 0
for(i in 1:nrow(viajes_22_sitp)){
#print(i)
det <- detailed_itineraries(r5r_core = r5r_core,
origins = viajes_22_origen[i,],
destinations = viajes_22_destino[i,],
mode = mode,
departure_datetime = departure_datetime,
max_walk_time = 180,
max_trip_duration = 300,
shortest_path = TRUE)
if(!is.na(det$total_distance[1])){
det <- det %>% filter(mode == "BUS")
viajes_22_sitp$dist_sitp[i] <- sum(det$distance)
}
}
print(Sys.time())
write.xlsx(viajes_22_sitp, "viajes_22_sitp.xlsx")
viajes_12_sitp <- read.xlsx("viajes_12_sitp.xlsx")
viajes_21_sitp <- read.xlsx("viajes_21_sitp.xlsx")
viajes_22_sitp <- read.xlsx("viajes_22_sitp.xlsx")
viajes_12 <- viajes_12 %>%
left_join(viajes_12_sitp[,c("cod_vj", "dist_sitp")], by = c("cod_vj"))
viajes_21 <- viajes_21 %>%
left_join(viajes_21_sitp[,c("cod_vj", "dist_sitp")], by = c("cod_vj"))
viajes_22 <- viajes_22 %>%
left_join(viajes_22_sitp[,c("cod_vj", "dist_sitp")], by = c("cod_vj"))
viajes_dist <- rbind(viajes_12, viajes_21, viajes_22)
viajes_dist$Distancia <- viajes_dist$dist_car
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "TransMilenio"] <- viajes_dist$dist_transmi[viajes_dist$modo_principal_typo == "TransMilenio"]
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "SITP Zonal"] <- viajes_dist$dist_sitp[viajes_dist$modo_principal_typo == "SITP Zonal"]
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "A pie"] <- viajes_dist$dist_walk[viajes_dist$modo_principal_typo == "A pie"]
viajes_dist$Distancia[viajes_dist$modo_principal_typo %in% c("Bicicleta")] <- viajes_dist$dist_bike[viajes_dist$modo_principal_typo %in% c("Bicicleta")]
#viajes_dist$Distancia <- as.numeric(viajes_dist$Distancia)
viajes_dist <- read.xlsx("viajes_dist_no_corr.xlsx")
Para los viajes en Transmilenio, dado que la distancia desde el punto de partida hasta la estación, y desde la estación hasta el punto de llegada, pueden ser largas, calculamos además la distancia por la red vial con ponderación peatón entre estos puntos, y la vamos sumando a la distancia calculada por la red troncal. El cálculo se realiza entonces con dodgr. Lo mismo sucede para los viajes con el SITP Zonal.
plot(net_transmi)
#Conversion de nudos de la red en formato sf. Incluye 140 estaciones + 25 nudos correspondientes a extremidades de ejes que no pudieron ser evitados en la construcción de la red
acceso_transmi <- net_transmi %>%
activate("nodes") %>%
st_geometry() %>%
st_as_sf()
#Matriz de distancias a pie entre hogares y estaciones
mat_walk_hog_transmi <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(acceso_transmi)), shortest = TRUE)
rownames(mat_walk_hog_transmi) <- Hogares$cod_hog
colnames(mat_walk_hog_transmi) <- rownames(acceso_transmi)
#Búsqueda de la estación de Transmilenio más cercana a cada hogar
distancia <- apply(mat_walk_hog_transmi[,1:dim(mat_walk_hog_transmi)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_hog_transmi)[apply(mat_walk_hog_transmi, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_hogar <- cbind(name, distancia)
estacion_mas_cerca_hogar$cod_hog <- as.numeric(rownames(estacion_mas_cerca_hogar))
#Matriz de distancias a pie entre centroides de ZAT y estaciones
mat_walk_centroid_transmi <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(acceso_transmi)), shortest = TRUE)
rownames(mat_walk_centroid_transmi) <- ZAT$ZAT
colnames(mat_walk_centroid_transmi) <- rownames(acceso_transmi)
#Búsqueda de la estación de Transmilenio más cercana a cada centroide de ZAT
distancia <- apply(mat_walk_centroid_transmi[,1:dim(mat_walk_centroid_transmi)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_centroid_transmi)[apply(mat_walk_centroid_transmi, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_centroid <- cbind(name, distancia)
estacion_mas_cerca_centroid$ZAT <- as.numeric(rownames(estacion_mas_cerca_centroid))
#Adición de la distancia hasta y desde Transmilenio en los viajes con este modo
Hog_TM <- viajes_dist %>% left_join(estacion_mas_cerca_hogar[,c("cod_hog", "distancia")], by = "cod_hog") %>%
group_by(cod_hog) %>%
summarize(dist_Hog_TM = unique(distancia), dist_TM_Hog = unique(distancia))
ZAT_TM <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_ori" = "ZAT")) %>%
group_by(zat_ori) %>%
summarize(dist_ZAT_TM = unique(distancia))
TM_ZAT <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_des" = "ZAT")) %>%
group_by(zat_des) %>%
summarize(dist_TM_ZAT = unique(distancia))
viajes_dist <- viajes_dist %>%
left_join(Hog_TM, by = "cod_hog") %>%
left_join(ZAT_TM, by = "zat_ori") %>%
left_join(TM_ZAT, by = "zat_des")
viajes_dist$dist_hasta_TM[viajes_dist$test_origen_destino == "HOGAR_Otro"] <- viajes_dist$dist_Hog_TM[viajes_dist$test_origen_destino == "HOGAR_Otro"]
viajes_dist$dist_hasta_TM[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")] <- viajes_dist$dist_ZAT_TM[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")]
viajes_dist$dist_desde_TM[viajes_dist$test_origen_destino == "Otro_HOGAR"] <- viajes_dist$dist_TM_Hog[viajes_dist$test_origen_destino == "Otro_HOGAR"]
viajes_dist$dist_desde_TM[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")] <- viajes_dist$dist_TM_ZAT[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")]
#Se trata de una solución provisional para calcular la distancia a la estación del SITP más cercana mientras resolvemos el problema con r5r.
#Paraderos oficiales del SITP. Incluye 8 517 paraderos del SITP
paraderos_SITP <- st_read(dsn = ".", layer = "Paraderos_SITP") %>% st_transform(4326)
## Reading layer `Paraderos_SITP' from data source
## `C:\Users\hugot\Documents\Rennes 2\Calcul distance SIG\Reprise perso\RMD Bogota 2023'
## using driver `ESRI Shapefile'
## Simple feature collection with 8517 features and 6 fields
## Geometry type: POINT
## Dimension: XY
## Bounding box: xmin: 587059 ymin: 486191.3 xmax: 610279.6 ymax: 532670.2
## Projected CRS: WGS 84 / UTM zone 18N
plot(st_geometry(paraderos_SITP))
#Matriz de distancias a pie entre hogares y paraderos
mat_walk_hog_SITP <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(Hogares)), to = as.data.frame(st_coordinates(paraderos_SITP)), shortest = TRUE)
rownames(mat_walk_hog_SITP) <- Hogares$cod_hog
colnames(mat_walk_hog_SITP) <- rownames(paraderos_SITP)
#Búsqueda del paradero de SITP más cercano a cada hogar
distancia <- apply(mat_walk_hog_SITP[,1:dim(mat_walk_hog_SITP)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_hog_SITP)[apply(mat_walk_hog_SITP, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_hogar <- cbind(name, distancia)
estacion_mas_cerca_hogar$cod_hog <- as.numeric(rownames(estacion_mas_cerca_hogar))
#Matriz de distancias a pie entre centroides de ZAT y paraderos
mat_walk_centroid_SITP <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATCentroids)), to = as.data.frame(st_coordinates(paraderos_SITP)), shortest = TRUE)
rownames(mat_walk_centroid_SITP) <- ZAT$ZAT
colnames(mat_walk_centroid_SITP) <- rownames(paraderos_SITP)
#Búsqueda del paradero de SITP más cercano a cada centroide de ZAT
distancia <- apply(mat_walk_centroid_SITP[,1:dim(mat_walk_centroid_SITP)[2]], FUN = min, MARGIN = 1)
name <- data.frame(name = colnames(mat_walk_centroid_SITP)[apply(mat_walk_centroid_SITP, 1, function(x) which(x == min(x))[1])])
estacion_mas_cerca_centroid <- cbind(name, distancia)
estacion_mas_cerca_centroid$ZAT <- as.numeric(rownames(estacion_mas_cerca_centroid))
#Adición de la distancia hasta y desde el SITP en los viajes con este modo
Hog_SITP <- viajes_dist %>% left_join(estacion_mas_cerca_hogar[,c("cod_hog", "distancia")], by = "cod_hog") %>%
group_by(cod_hog) %>%
summarize(dist_Hog_SITP = unique(distancia), dist_SITP_Hog = unique(distancia))
ZAT_SITP <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_ori" = "ZAT")) %>%
group_by(zat_ori) %>%
summarize(dist_ZAT_SITP = unique(distancia))
SITP_ZAT <- viajes_dist %>% left_join(estacion_mas_cerca_centroid[,c("ZAT", "distancia")], by = c("zat_des" = "ZAT")) %>%
group_by(zat_des) %>%
summarize(dist_SITP_ZAT = unique(distancia))
viajes_dist <- viajes_dist %>%
left_join(Hog_SITP, by = "cod_hog") %>%
left_join(ZAT_SITP, by = "zat_ori") %>%
left_join(SITP_ZAT, by = "zat_des")
viajes_dist$dist_hasta_SITP[viajes_dist$test_origen_destino == "HOGAR_Otro"] <- viajes_dist$dist_Hog_SITP[viajes_dist$test_origen_destino == "HOGAR_Otro"]
viajes_dist$dist_hasta_SITP[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")] <- viajes_dist$dist_ZAT_SITP[viajes_dist$test_origen_destino %in% c("Otro_HOGAR", "Otro_Otro")]
viajes_dist$dist_desde_SITP[viajes_dist$test_origen_destino == "Otro_HOGAR"] <- viajes_dist$dist_SITP_Hog[viajes_dist$test_origen_destino == "Otro_HOGAR"]
viajes_dist$dist_desde_SITP[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")] <- viajes_dist$dist_SITP_ZAT[viajes_dist$test_origen_destino %in% c("HOGAR_Otro", "Otro_Otro")]
Hallamos la distancia promedio desde cada red hasta los puntos de la matriz (con resolución 200m).
comparacion <- data.frame(red = c("Transmilenio - Hogar","Transmilenio - ZAT","SITP - Hogar", "SITP - ZAT"),
dist_promedio_desde_matriz = c(mean(Hog_TM$dist_Hog_TM),
mean(ZAT_TM$dist_ZAT_TM),
mean(Hog_SITP$dist_Hog_SITP),
mean(ZAT_SITP$dist_ZAT_SITP)),
dist_mediana_desde_matriz = c(median(Hog_TM$dist_Hog_TM),
median(ZAT_TM$dist_ZAT_TM),
median(Hog_SITP$dist_Hog_SITP),
median(ZAT_SITP$dist_ZAT_SITP)))
colnames(comparacion) <- c("Red", "Promedio", "Mediana")
kable(comparacion, caption = "Comparación de las distancias entre los Hogares, los centroides de ZAT y las redes de transporte colectivo", format.args = list(big.mark = ","))
| Red | Promedio | Mediana |
|---|---|---|
| Transmilenio - Hogar | 5,671.487 | 1,941.5376 |
| Transmilenio - ZAT | 5,408.001 | 1,912.4224 |
| SITP - Hogar | 2,917.074 | 213.0079 |
| SITP - ZAT | 2,779.874 | 249.1494 |
Para los trayectos intrazonas, dos metodos fueron probados en la base de la EODH 2019:
En nuestro caso, escogemos el método “de la media” que tiende a dar estimaciones superiores al método “del CEREMA”.
Para el caso de los viajes intrazonas en Transmilenio, el criterio para detectarlos es que la distancia por la red de Transmilenio sea nula.
#Seleccion de las ZAT que tienen trayectos internos
viajes_intra <- viajes_dist[viajes_dist$zat_ori == viajes_dist$zat_des & viajes_dist$Distancia == 0,c("zat_ori", "zat_des", "cod_vj", "modo_principal_typo")]
viajes_intra$id_row <- 1:nrow(viajes_intra)
ZAT_intra <- ZAT[ZAT$ZAT %in% as.numeric(unique(viajes_intra$zat_des)),]
#Tomando muestra en cada ZAT seleccionada
#Tamaño de la muestra
s <- 10
d <- nrow(ZAT_intra)
set.seed(1)
ZATSampled <- st_sample(ZAT_intra, rep(s,d))
ZATSampled <- st_sf(ZATSampled) %>%
st_join(ZAT[,"ZAT"],
join = st_intersects)
ZATSampled$row <- rep(c(1:s),d)
#Calculando la matriz de distancias por la red vial, ponderación peatón
mat_intra <- dodgr_dists(graph = net_walk, from = as.data.frame(st_coordinates(ZATSampled)), to = as.data.frame(st_coordinates(ZATSampled)), shortest = TRUE)
rownames(mat_intra) <- ZATSampled$ZAT
colnames(mat_intra) <- ZATSampled$ZAT
comb <- function(z) {
merge(viajes_intra, z, by = c("zat_ori","zat_des")) %>%
arrange(id_row) %>%
dplyr::select(c(-zat_ori,-zat_des, -cod_vj, -modo_principal_typo, -id_row))
}
doFuture::registerDoFuture()
future::plan(multisession, workers = 8)
options(future.globals.maxSize= 1048576000)
viajes2 <- foreach(u = 1:s, .combine = 'cbind') %:%
foreach(v = 1:s, .combine = 'cbind') %dofuture% {
subset_i = seq(from = u, to = s*d, by = s)
subset_j = seq(from = v, to = s*d, by = s)
z <- mat_intra[subset_i,subset_j]
colnames(z) <- ZAT_intra$ZAT
rownames(z) <- ZAT_intra$ZAT
z <- data.frame(rownames(as.matrix(z))[row(as.matrix(z))],
colnames(as.matrix(z))[col(as.matrix(z))],
c(as.matrix(z)))
colnames(z) <- c("zat_ori", "zat_des", paste0("dist",(u-1)*s+(v-1)))
x <- comb(z)
}
viajes_intra <- cbind(viajes_intra,viajes2)
#índices de las columnas con distancia interna igual a cero
index_zero <- rep("VIDE",s)
for(i in 1:s){
index_zero[i] <- paste0("dist",(i-1)*s+(i-1))
}
#Distancia interna = promedio de las distancias excluyendo los ceros.
viajes_intra_non_zero <- viajes_intra[,!(colnames(viajes_intra) %in% index_zero)]
viajes_intra$dist_mean_non_zeros <- rowMeans(viajes_intra_non_zero[,6:95])
#Sustituyendo los ceros en la base de distancias
viajes_dist$Distancia[viajes_dist$zat_ori == viajes_dist$zat_des & viajes_dist$Distancia == 0] <- viajes_intra$dist_mean_non_zeros
Identificamos los viajes de ida y vuelta. Un viaje de ida y vuelta se identifica porque cumple las siguientes características:
Nota: dado que no pudimos identificar todos los viajes que inician en el hogar, sino solo el primero del día, este algoritmo solo puede detectar máximo una ida y vuelta por persona
#Esta etapa es necesaria en el marco del control de velocidad de la siguiente etapa. Se define como ida y vuelta un trayecto casa -> cierto lugar -> casa realizado con el mismo modo de transporte. El lugar de destino del primer viaje y el lugar de origen del segundo deben coincidir. Si el modo de transporte cambia, no lo definimos como ida y vuelta.
ida <- viajes_dist[viajes_dist$lugar_origen == "Hogar", c("cod_pers", "cod_vj", "lugar_destino", "modo_principal_typo")]
vuelta <- viajes_dist[viajes_dist$lugar_destino == "Hogar", c("cod_pers", "cod_vj", "lugar_origen", "modo_principal_typo")]
ida_vuelta <- na.omit(ida %>% left_join(vuelta, by = "cod_pers"))
ida_vuelta$is_ida_vuelta <- (ida_vuelta$lugar_destino == ida_vuelta$lugar_origen) & (ida_vuelta$modo_principal_typo.x == ida_vuelta$modo_principal_typo.y)
ida <- ida_vuelta[,c("cod_vj.x", "is_ida_vuelta")]
colnames(ida)[colnames(ida) == "cod_vj.x"] <- "cod_vj"
ida <- ida %>% group_by(cod_vj) %>% summarize(is_ida_vuelta = as.logical(max(is_ida_vuelta)))
vuelta <- ida_vuelta[,c("cod_vj.y", "is_ida_vuelta")]
colnames(vuelta)[colnames(vuelta) == "cod_vj.y"] <- "cod_vj"
vuelta <- vuelta %>% group_by(cod_vj) %>% summarize(is_ida_vuelta = as.logical(max(is_ida_vuelta)))
viajes_dist <- viajes_dist %>% left_join(ida, by = c("cod_vj"))
viajes_dist <- viajes_dist %>% left_join(vuelta, by = c("cod_vj"))
viajes_dist$is_ida_vuelta.x[is.na(viajes_dist$is_ida_vuelta.x)] <- FALSE
viajes_dist$is_ida_vuelta.y[is.na(viajes_dist$is_ida_vuelta.y)] <- FALSE
viajes_dist$is_ida_vuelta <- viajes_dist$is_ida_vuelta.x | viajes_dist$is_ida_vuelta.y
viajes_dist$is_ida_vuelta.x <- NULL
viajes_dist$is_ida_vuelta.y <- NULL
#velocidad en km/h
viajes_dist$velocidad <- 60*viajes_dist$Distancia/(1000*viajes_dist$duracion_min)
#8139 viajes a pie son demasiado rápidos (23%)
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "A pie" & viajes_dist$velocidad > 10] <- 1000*viajes_dist$duracion_min[viajes_dist$modo_principal_typo == "A pie" & viajes_dist$velocidad > 10]*10/60
#250 viajes en bici son demasiado rápidos (5%)
viajes_dist$Distancia[viajes_dist$modo_principal_typo == "Bicicleta" & viajes_dist$velocidad > 30] <- 1000*viajes_dist$duracion_min[viajes_dist$modo_principal_typo == "Bicicleta" & viajes_dist$velocidad > 30]*30/60
#963 viajes en modos motorizados son demasiado rápidos (2%)
viajes_dist$Distancia[!(viajes_dist$modo_principal_typo %in% c("A pie","Bicicleta")) & viajes_dist$velocidad > 50] <- 1000*viajes_dist$duracion_min[!(viajes_dist$modo_principal_typo %in% c("A pie","Bicicleta")) & viajes_dist$velocidad > 50]*50/60
#se aplica una ultima correccion: los viajes de ida y vuelta deben tener la misma distancia. Se corrige la distancia de los viajes demasiado rápidos, pero en el caso de los viajes de ida y vuelta, dado que esperamos que la distancia coincida en la ida y el regreso, le vamos a asignar el maximo entre los dos valores.
ida_vuelta <- viajes_dist[viajes_dist$is_ida_vuelta == TRUE,] %>%
group_by(cod_pers, modo_principal_typo) %>%
summarise(Distancia_AR = max(Distancia))
viajes_dist <- viajes_dist %>% left_join(ida_vuelta, by = c("cod_pers", "modo_principal_typo"))
viajes_dist$Distancia[viajes_dist$is_ida_vuelta == TRUE] <- viajes_dist$Distancia_AR[viajes_dist$is_ida_vuelta == TRUE]
viajes_dist$Distancia_AR <- NULL
viajes_dist$velocidad <- NULL
viajes_dist$Distancia_incl_caminata <- viajes_dist$Distancia
viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo == "TransMilenio"] <- viajes_dist$Distancia[viajes_dist$modo_principal_typo == "TransMilenio"] + viajes_dist$dist_hasta_TM[viajes_dist$modo_principal_typo == "TransMilenio"] + viajes_dist$dist_desde_TM[viajes_dist$modo_principal_typo == "TransMilenio"]
viajes_dist$Distancia_incl_caminata[viajes_dist$modo_principal_typo == "SITP Zonal"] <- viajes_dist$Distancia[viajes_dist$modo_principal_typo == "SITP Zonal"] + viajes_dist$dist_hasta_SITP[viajes_dist$modo_principal_typo == "SITP Zonal"] + viajes_dist$dist_desde_SITP[viajes_dist$modo_principal_typo == "SITP Zonal"]
Por fin podemos crear la base final con las variables correspondientes a distancia, emisiones de gases de efecto invernadero y contaminantes aéreos.
viajes_dist <- viajes_dist %>% left_join(EF, by = c("modo_principal_typo" = "modo_principal"))
viajes_dist$`CO2-eq` <- viajes_dist$`CO2-eq`*viajes_dist$Distancia/1000
viajes_dist$CO <- viajes_dist$CO*viajes_dist$Distancia/1000
viajes_dist$NOx <- viajes_dist$NOx*viajes_dist$Distancia/1000
viajes_dist$SO2 <- viajes_dist$SO2*viajes_dist$Distancia/1000
viajes_dist$COV <- viajes_dist$COV*viajes_dist$Distancia/1000
viajes_dist$PM.2.5 <- viajes_dist$PM.2.5*viajes_dist$Distancia/1000
viajes_dist$PM.10 <- viajes_dist$PM.10*viajes_dist$Distancia/1000
#write.xlsx(viajes_dist, paste0("viajes_dist_final_con_caminata_hasta_TM.xlsx"), overwrite = TRUE)
viajes_dist <- read.xlsx("viajes_dist_final_con_caminata_hasta_TM.xlsx")
estadisticas <- viajes_dist %>%
group_by(modo_principal_typo) %>%
summarize(numero_viajes = round(sum(fexp_vj),0),
mpkt = round(sum(fexp_vj*Distancia/1000000000),1),
mpkt_incl_caminata = round(sum(fexp_vj*Distancia_incl_caminata/1000000000),1),
reparto_modal = round(100*sum(fexp_vj)/sum(viajes_dist$fexp_vj),1),
reparto_modal_pkt = round(100*sum(fexp_vj*Distancia_incl_caminata)/sum(viajes_dist$fexp_vj*viajes_dist$Distancia_incl_caminata),1),
`CO2-eq`= round(sum(fexp_vj*`CO2-eq`/1000000),1),
dist_por_viaj = round(sum(fexp_vj*Distancia_incl_caminata/1000)/sum(fexp_vj),1))
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]
total <- data.frame("TOTAL",
sum(estadisticas$numero_viajes),
sum(estadisticas$mpkt),
sum(estadisticas$mpkt_incl_caminata),
sum(estadisticas$reparto_modal),
sum(estadisticas$reparto_modal_pkt),
sum(estadisticas$`CO2-eq`),
round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))
colnames(total) <- colnames(estadisticas)
estadisticas <- rbind(estadisticas, total)
colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")
kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))
| Modo principal | Número de viajes | Millón PKT (solo tramo principal) | Millón PKT (total incluyendo caminata inicial y final) | % Número de viajes | % PKT | Emisiones GEI (tCO2-eq) | Distancia media por viaje (km) |
|---|---|---|---|---|---|---|---|
| A pie | 6,134,682 | 10.8 | 10.8 | 37.4 | 9.1 | 0.0 | 1.8 |
| SITP Zonal | 2,237,192 | 24.6 | 25.9 | 13.6 | 22.0 | 1,469.7 | 11.6 |
| TransMilenio | 2,008,743 | 26.3 | 34.0 | 12.3 | 28.8 | 105.8 | 16.9 |
| Auto | 1,916,020 | 16.0 | 16.0 | 11.7 | 13.6 | 2,433.2 | 8.4 |
| Moto | 1,053,826 | 10.9 | 10.9 | 6.4 | 9.3 | 834.5 | 10.4 |
| Bicicleta | 1,037,486 | 5.7 | 5.7 | 6.3 | 4.9 | 0.0 | 5.5 |
| Transporte publico individual | 670,709 | 4.3 | 4.3 | 4.1 | 3.6 | 695.8 | 6.4 |
| Intermunicipal | 432,446 | 4.1 | 4.1 | 2.6 | 3.5 | 293.7 | 9.6 |
| Transporte Escolar | 293,703 | 1.9 | 1.9 | 1.8 | 1.6 | 96.4 | 6.4 |
| Otro | 250,446 | 2.1 | 2.1 | 1.5 | 1.8 | 106.0 | 8.3 |
| Alimentador | 178,527 | 1.1 | 1.1 | 1.1 | 0.9 | 9.2 | 6.1 |
| Transporte informal | 119,618 | 0.8 | 0.8 | 0.7 | 0.7 | 64.0 | 7.0 |
| Bicitaxi | 47,601 | 0.2 | 0.2 | 0.3 | 0.2 | 10.5 | 4.8 |
| Patineta | 12,677 | 0.1 | 0.1 | 0.1 | 0.1 | 0.2 | 5.9 |
| Cable | 2,269 | 0.0 | 0.0 | 0.0 | 0.0 | 0.1 | 3.5 |
| TOTAL | 16,395,945 | 108.9 | 117.9 | 99.9 | 100.1 | 6,119.1 | 7.2 |
#write.xlsx(estadisticas, "resumen_distancias_emisiones_con_caminata_hasta_TM.xlsx", overwrite = TRUE)
estadisticas <- viajes_dist %>%
group_by(modo_principal_comparado) %>%
summarize(numero_viajes = round(sum(fexp_vj),0),
mpkt = round(sum(fexp_vj*Distancia/1000000000),1),
mpkt_incl_caminata = round(sum(fexp_vj*Distancia_incl_caminata/1000000000),1),
reparto_modal = round(100*sum(fexp_vj)/sum(viajes_dist$fexp_vj),1),
reparto_modal_pkt = round(100*sum(fexp_vj*Distancia_incl_caminata)/sum(viajes_dist$fexp_vj*viajes_dist$Distancia_incl_caminata),1),
`CO2-eq`= round(sum(fexp_vj*`CO2-eq`/1000000),1),
dist_por_viaj = round(sum(fexp_vj*Distancia_incl_caminata/1000)/sum(fexp_vj),1))
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]
total <- data.frame("TOTAL",
sum(estadisticas$numero_viajes),
sum(estadisticas$mpkt),
sum(estadisticas$mpkt_incl_caminata),
sum(estadisticas$reparto_modal),
sum(estadisticas$reparto_modal_pkt),
sum(estadisticas$`CO2-eq`),
round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))
colnames(total) <- colnames(estadisticas)
estadisticas <- rbind(estadisticas, total)
colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")
kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))
| Modo principal | Número de viajes | Millón PKT (solo tramo principal) | Millón PKT (total incluyendo caminata inicial y final) | % Número de viajes | % PKT | Emisiones GEI (tCO2-eq) | Distancia media por viaje (km) |
|---|---|---|---|---|---|---|---|
| A pie | 6,134,682 | 10.8 | 10.8 | 37.4 | 9.1 | 0.0 | 1.8 |
| Transporte publico | 4,859,177 | 56.1 | 65.1 | 29.6 | 55.2 | 1,878.5 | 13.4 |
| Auto | 1,916,020 | 16.0 | 16.0 | 11.7 | 13.6 | 2,433.2 | 8.4 |
| Moto | 1,053,826 | 10.9 | 10.9 | 6.4 | 9.3 | 834.5 | 10.4 |
| Bicicleta | 1,037,486 | 5.7 | 5.7 | 6.3 | 4.9 | 0.0 | 5.5 |
| Taxi / Carro por aplicación | 670,709 | 4.3 | 4.3 | 4.1 | 3.6 | 695.8 | 6.4 |
| Transporte Escolar | 293,703 | 1.9 | 1.9 | 1.8 | 1.6 | 96.4 | 6.4 |
| Otro | 263,123 | 2.1 | 2.1 | 1.6 | 1.8 | 106.1 | 8.2 |
| Transporte informal | 167,219 | 1.1 | 1.1 | 1.0 | 0.9 | 74.5 | 6.4 |
| TOTAL | 16,395,945 | 108.9 | 117.9 | 99.9 | 100.0 | 6,119.0 | 7.2 |
#write.xlsx(estadisticas, "comparado_distancias_emisiones_con_caminata_hasta_TM.xlsx", overwrite = TRUE)
estadisticas <- viajes_dist %>%
group_by(modo_principal_LATS) %>%
summarize(numero_viajes = round(sum(fexp_vj),0),
mpkt = round(sum(fexp_vj*Distancia/1000000000),1),
mpkt_incl_caminata = round(sum(fexp_vj*Distancia_incl_caminata/1000000000),1),
reparto_modal = round(100*sum(fexp_vj)/sum(viajes_dist$fexp_vj),1),
reparto_modal_pkt = round(100*sum(fexp_vj*Distancia_incl_caminata)/sum(viajes_dist$fexp_vj*viajes_dist$Distancia_incl_caminata),1),
`CO2-eq`= round(sum(fexp_vj*`CO2-eq`/1000000),1),
dist_por_viaj = round(sum(fexp_vj*Distancia_incl_caminata/1000)/sum(fexp_vj),1))
estadisticas <- estadisticas[order(estadisticas$numero_viajes, decreasing = TRUE),]
total <- data.frame("TOTAL",
sum(estadisticas$numero_viajes),
sum(estadisticas$mpkt),
sum(estadisticas$mpkt_incl_caminata),
sum(estadisticas$reparto_modal),
sum(estadisticas$reparto_modal_pkt),
sum(estadisticas$`CO2-eq`),
round(sum(estadisticas$numero_viajes*estadisticas$dist_por_viaj)/sum(estadisticas$numero_viajes),1))
colnames(total) <- colnames(estadisticas)
estadisticas <- rbind(estadisticas, total)
colnames(estadisticas) <- c("Modo principal","Número de viajes","Millón PKT (solo tramo principal)", "Millón PKT (total incluyendo caminata inicial y final)","% Número de viajes","% PKT","Emisiones GEI (tCO2-eq)", "Distancia media por viaje (km)")
kable(estadisticas, caption = "Estadísticas básicas", format.args = list(big.mark = ","))
| Modo principal | Número de viajes | Millón PKT (solo tramo principal) | Millón PKT (total incluyendo caminata inicial y final) | % Número de viajes | % PKT | Emisiones GEI (tCO2-eq) | Distancia media por viaje (km) |
|---|---|---|---|---|---|---|---|
| Walking | 6,134,682 | 10.8 | 10.8 | 37.4 | 9.1 | 0.0 | 1.8 |
| Regular bus | 2,751,356 | 28.6 | 29.9 | 16.8 | 25.4 | 1,704.2 | 10.9 |
| BRT | 2,107,820 | 27.5 | 35.2 | 12.9 | 29.8 | 174.3 | 16.7 |
| Private car | 1,916,020 | 16.0 | 16.0 | 11.7 | 13.6 | 2,433.2 | 8.4 |
| Moto | 1,053,826 | 10.9 | 10.9 | 6.4 | 9.3 | 834.5 | 10.4 |
| Bike | 1,037,486 | 5.7 | 5.7 | 6.3 | 4.9 | 0.0 | 5.5 |
| Taxi | 670,709 | 4.3 | 4.3 | 4.1 | 3.6 | 695.8 | 6.4 |
| Transporte Escolar | 293,703 | 1.9 | 1.9 | 1.8 | 1.6 | 96.4 | 6.4 |
| Other | 263,123 | 2.1 | 2.1 | 1.6 | 1.8 | 106.1 | 8.2 |
| Paratransit | 167,219 | 1.1 | 1.1 | 1.0 | 0.9 | 74.5 | 6.4 |
| TOTAL | 16,395,944 | 108.9 | 117.9 | 100.0 | 100.0 | 6,119.0 | 7.2 |
#write.xlsx(estadisticas, "LATS_distancias_emisiones_con_caminata_hasta_TM.xlsx", overwrite = TRUE)
Finalmente, producimos mapas por zonas (UTAM)
centroides_limites_bogota <- st_centroid(limites_bogota) %>%
mutate(lon = st_coordinates(.)[,1],
lat = st_coordinates(.)[,2])
centroides_limites_bogota$MUNI <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$MUNI[centroides_limites_bogota$NOMBR_MUNI == "BOGOTA"] <- ""
centroides_limites_bogota$LOCA <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$LOCA[centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"] <- ""
centroides_limites_bogota$DISPLAY <- centroides_limites_bogota$NOMBRE
centroides_limites_bogota$DISPLAY[!((centroides_limites_bogota$NOMBR_MUNI == "BOGOTA" & centroides_limites_bogota$NOMBRE %in% c("BOSA", "SANTA FE", "CIUDAD BOLIVAR", "TEUSAQUILLO", "FONTIBON", "ENGATIVA", "SUBA", "USAQUEN", "USME", "SAN CRISTOBAL", "CHAPINERO", "KENNEDY")) | (centroides_limites_bogota$NOMBR_MUNI != "BOGOTA"))] <- ""
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "SUBA"] <- 4.76
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "SOACHA"] <- 4.57
centroides_limites_bogota$lon[centroides_limites_bogota$NOMBRE == "CHAPINERO"] <- -74.027
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "USAQUEN"] <- 4.74
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "CHAPINERO"] <- 4.68
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "TEUSAQUILLO"] <- 4.65
centroides_limites_bogota$lat[centroides_limites_bogota$NOMBRE == "KENNEDY"] <- 4.635
centroides_limites_bogota$DISPLAY[centroides_limites_bogota$NOMBRE == "CIUDAD BOLIVAR"] <- "CIUDAD \nBOLIVAR"
ZAT_Poblacion <- Personas %>%
group_by(zat_hg) %>%
summarize(pop = sum(`fexp_per>5años`))
ZAT <- ZAT %>% left_join(ZAT_Poblacion, by = c("ZAT" = "zat_hg"))
ZAT$pop[is.na(ZAT$pop)] <- 0
UTAM <- UTAM %>% st_transform(4326)
Dist_ZAT <- viajes_dist[,c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
group_by(zat_hg) %>%
summarize(dist = sum(fexp_vj*Distancia_incl_caminata/1000))
ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))
ZAT$dist[is.na(ZAT$dist)] <- 0
UTAM_indic <- ZAT %>%
st_drop_geometry() %>%
group_by(UTAM) %>%
summarise(dist = sum(dist),
pop = sum(pop),
dist_capita = dist/pop)
UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")
#Colores
bks_total <- round(as.numeric(quantile(UTAM$dist_capita, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_total,max(UTAM$dist_capita))
UTAM$binned_dist_capita <- cut(UTAM$dist_capita,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_dist_capita))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(6,"Purples")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Distancia cotidiana promedio per \ncápita todos modos incluidos (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
Dist_ZAT <- viajes_dist[viajes_dist$modo_principal_comparado == "A pie",c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
group_by(zat_hg) %>%
summarize(dist_walk = sum(fexp_vj*Distancia_incl_caminata/1000))
ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))
ZAT$dist_walk[is.na(ZAT$dist_walk)] <- 0
UTAM_indic <- ZAT %>%
st_drop_geometry() %>%
group_by(UTAM) %>%
summarise(dist_walk = sum(dist_walk),
walk_intensity = sum(dist_walk)/sum(pop))
UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")
#Colores
bks_walk <- round(as.numeric(quantile(UTAM$walk_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_walk,max(UTAM$walk_intensity))
UTAM$binned_walk_intensity <- cut(UTAM$walk_intensity,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_walk_intensity))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(6,"Reds")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Distancia cotidiana promedio \nper cápita caminando (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
Dist_ZAT <- viajes_dist[viajes_dist$modo_principal_comparado == "Auto",c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
group_by(zat_hg) %>%
summarize(dist_auto = sum(fexp_vj*Distancia_incl_caminata/1000))
ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))
ZAT$dist_auto[is.na(ZAT$dist_auto)] <- 0
UTAM_indic <- ZAT %>%
st_drop_geometry() %>%
group_by(UTAM) %>%
summarise(dist_auto = sum(dist_auto),
auto_intensity = sum(dist_auto)/sum(pop))
UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")
bks_car <- round(as.numeric(quantile(UTAM$auto_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),2)
bks <- c(bks_car,max(UTAM$auto_intensity))
UTAM$binned_auto_intensity <- cut(UTAM$auto_intensity,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_auto_intensity))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(6,"RdPu")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Distancia cotidiana promedio per \ncápita en auto particular (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
Dist_ZAT <- viajes_dist[viajes_dist$modo_principal_comparado == "Transporte publico",c("zat_hg", "modo_principal_comparado", "Distancia_incl_caminata", "fexp_vj")]%>%
left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
group_by(zat_hg) %>%
summarize(dist_transit = sum(fexp_vj*Distancia_incl_caminata/1000))
ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))
ZAT$dist_transit[is.na(ZAT$dist_transit)] <- 0
UTAM_indic <- ZAT %>%
st_drop_geometry() %>%
group_by(UTAM) %>%
summarise(dist_transit = sum(dist_transit),
transit_intensity = sum(dist_transit)/sum(pop))
UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")
bks_transit <- round(as.numeric(quantile(UTAM$transit_intensity, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),3)
bks <- c(bks_transit,max(UTAM$transit_intensity))
UTAM$binned_transit_intensity <- cut(UTAM$transit_intensity,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_transit_intensity))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(6,"Greens")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Distancia cotidiana promedio per \ncápita en transporte público (km)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
Dist_ZAT <- viajes_dist[,c("zat_hg", "modo_principal_comparado", "CO2-eq", "fexp_vj")]%>%
left_join(ZAT[,c("ZAT")], by = c("zat_hg" = "ZAT")) %>%
group_by(zat_hg) %>%
summarize(co2 = sum(fexp_vj*`CO2-eq`/1000000))
ZAT <- ZAT %>% left_join(Dist_ZAT, by = c("ZAT" = "zat_hg"))
ZAT$co2[is.na(ZAT$co2)] <- 0
UTAM_indic <- ZAT %>%
st_drop_geometry() %>%
group_by(UTAM) %>%
summarise(co2 = sum(co2),
co2_capita = 1000*sum(co2)/sum(pop))
UTAM <- UTAM %>% left_join(UTAM_indic, by = "UTAM")
bks_co2 <- round(as.numeric(quantile(UTAM$co2_capita, probs = c(0,0.1,0.25,0.5,0.75,0.9), na.rm = TRUE)),5)
bks <- c(bks_co2,max(UTAM$co2_capita))
UTAM$binned_co2_capita <- cut(UTAM$co2_capita,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_co2_capita))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Límite de municipio fuera de Bogotá DC \nLímite de localidad dentro de Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(6,"Oranges")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Emisiones cotidianas promedio \nper cápita (kg CO2-eq)", caption = "Autor: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nFuente: Encuesta Origen-Destino de Hogares, Secretaría Distrital de Movilidad, 2023 \nDiscretización con los cuantiles 0.10, 0.25, 0.50, 0.75, 0.90", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 9, face = "italic", hjust = 0, vjust = 90),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
#Colores
bks_total <- round(as.numeric(quantile(UTAM$dist_capita, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_total,max(UTAM$dist_capita))
UTAM$binned_dist_capita <- cut(UTAM$dist_capita,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_dist_capita))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 4.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(5,"Purples")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Average daily distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
#Colores
bks_walk <- round(as.numeric(quantile(UTAM$walk_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_walk,max(UTAM$walk_intensity))
UTAM$binned_walk_intensity <- cut(UTAM$walk_intensity,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_walk_intensity))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(5,"Reds")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Average daily walking distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
bks_car <- round(as.numeric(quantile(UTAM$auto_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),2)
bks <- c(bks_car,max(UTAM$auto_intensity))
UTAM$binned_auto_intensity <- cut(UTAM$auto_intensity,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_auto_intensity))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(5,"RdPu")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Average daily private car distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
bks_transit <- round(as.numeric(quantile(UTAM$transit_intensity, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),3)
bks <- c(bks_transit,max(UTAM$transit_intensity))
UTAM$binned_transit_intensity <- cut(UTAM$transit_intensity,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_transit_intensity))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(5,"Greens")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Average daily public transport distance \nper capita (km)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
bks_co2 <- round(as.numeric(quantile(UTAM$co2_capita, probs = c(0,0.2,0.4,0.6,0.8), na.rm = TRUE)),5)
bks <- c(bks_co2,max(UTAM$co2_capita))
UTAM$binned_co2_capita <- cut(UTAM$co2_capita,bks)
ggplot()+
theme_bw()+
geom_sf(data = UTAM, color = NA, aes(fill = binned_co2_capita))+
geom_sf(data = limites_bogota, fill = NA, aes(linetype = "Municipality boundary outside Bogotá DC \nLocality boundary inside Bogotá DC"), color = "grey")+
#geom_text(data= centroides_limites_bogota,aes(x=lon, y=lat, label = DISPLAY), color = "black", fontface = "plain", size = 2.5, check_overlap = FALSE) +
scale_fill_manual(values = brewer.pal(5,"Oranges")) +
coord_sf(xlim = lon_bounds_Bogota, ylim = lat_bounds_Bogota, datum = NA)+
labs(title = "", fill = "Average daily GHG emissions \nper capita (kg CO2-eq)", caption = "Author: Hugo Thomas - Université Rennes 2 / Universidad de los Andes \nSource: Households Mobility Survey, Bogotá Mobility Department (SDM), 2023 \nDiscretization in quintiles", linetype = "") +
theme(legend.position = "bottom",
legend.direction = "vertical",
plot.title = element_text(size = 15, face = "bold", hjust = 0.5, vjust = 0),
legend.title = element_text(size = 13),
legend.text = element_text(size=13),
plot.caption = element_text(size = 11, face = "italic", hjust = 0, vjust = 70),
legend.text.align = 1)+
guides(linetype = guide_legend(label.hjust = 0, order = 1))+
labs(x = "", y = "") +
annotation_scale(location = "br", height = unit(0.2, "cm"),text_cex = 0.8) +
annotation_north_arrow(location = "tl", which_north = "true", height = unit(1, "cm"), width = unit(1, "cm"), style = north_arrow_orienteering(text_size = 8))
setwd("C:/Users/hugot/Documents/Uniandes/LATS")
UTAM_ICS <- read.xlsx("UTAM_ICS.xlsx")
UTAM <- UTAM %>% left_join(UTAM_ICS[, c("UTAM", "propICS12")], by = "UTAM")
correlaciones <- rbind(data.frame(cor.test(UTAM$propICS12, UTAM$auto_intensity, method = "pearson")[c(4,3)]),
data.frame(cor.test(UTAM$walk_intensity, UTAM$auto_intensity, method = "pearson")[c(4,3)]),
data.frame(cor.test(UTAM$walk_intensity, UTAM$transit_intensity, method = "pearson")[c(4,3)]),
data.frame(cor.test(UTAM$propICS12, UTAM$walk_intensity, method = "pearson")[c(4,3)]),
data.frame(cor.test(UTAM$propICS12, UTAM$transit_intensity, method = "pearson")[c(4,3)]),
data.frame(cor.test(UTAM$co2_capita, UTAM$auto_intensity, method = "pearson")[c(4,3)]),
data.frame(cor.test(UTAM$propICS12, UTAM$co2_capita, method = "pearson")[c(4,3)]))
rownames(correlaciones) <- c("auto_ICS",
"walk_auto",
"walk_transit",
"ICS_walk",
"ICS_transit",
"GES_auto",
"ICS_GES")
kable(correlaciones, caption = "Correlaciones", format.args = list(big.mark = ","))
| estimate | p.value | |
|---|---|---|
| auto_ICS | -0.4873969 | 0.0000000 |
| walk_auto | -0.4482729 | 0.0000001 |
| walk_transit | 0.4210982 | 0.0000005 |
| ICS_walk | 0.4398057 | 0.0000001 |
| ICS_transit | 0.5627718 | 0.0000000 |
| GES_auto | 0.8709694 | 0.0000000 |
| ICS_GES | -0.2907798 | 0.0007185 |
correlaciones <- cbind(rownames(correlaciones), correlaciones)
write.xlsx(correlaciones, "correlaciones.xlsx", overwrite = TRUE)
UTAM_top_10 <- UTAM[order(UTAM$propICS12, decreasing = FALSE),]
select <-as.integer(0.1*nrow(UTAM_top_10))
UTAM_top_10 <- UTAM_top_10[1:select,]
UTAM_bottom_40 <- UTAM[order(UTAM$propICS12, decreasing = TRUE),]
select <-as.integer(0.4*nrow(UTAM_bottom_40))
UTAM_bottom_40 <- UTAM_bottom_40[1:select,]
palma <- rbind(data.frame(top_10 = weighted.mean(x = UTAM_top_10$dist_capita, w = UTAM_top_10$pop),
bottom_40 = weighted.mean(x = UTAM_bottom_40$dist_capita, w = UTAM_bottom_40$pop)),
data.frame(top_10 = weighted.mean(x = UTAM_top_10$auto_intensity, w = UTAM_top_10$pop),
bottom_40 = weighted.mean(x = UTAM_bottom_40$auto_intensity, w = UTAM_bottom_40$pop)),
data.frame(top_10 = weighted.mean(x = UTAM_top_10$walk_intensity, w = UTAM_top_10$pop),
bottom_40 = weighted.mean(x = UTAM_bottom_40$walk_intensity, w = UTAM_bottom_40$pop)),
data.frame(top_10 = weighted.mean(x = UTAM_top_10$transit_intensity, w = UTAM_top_10$pop),
bottom_40 = weighted.mean(x = UTAM_bottom_40$transit_intensity, w = UTAM_bottom_40$pop)),
data.frame(top_10 = weighted.mean(x = UTAM_top_10$co2_capita, w = UTAM_top_10$pop),
bottom_40 = weighted.mean(x = UTAM_bottom_40$co2_capita, w = UTAM_bottom_40$pop)))
palma$palma <- palma$top_10/palma$bottom_40
rownames(palma) <- c("dist_capita",
"auto_intensity",
"walk_intensity",
"transit_intensity",
"co2_capita")
kable(palma, caption = "Palma ratios", format.args = list(big.mark = ","))
| top_10 | bottom_40 | palma | |
|---|---|---|---|
| dist_capita | 10.4797194 | 14.1213597 | 0.7421183 |
| auto_intensity | 3.6798541 | 1.0520899 | 3.4976612 |
| walk_intensity | 0.7930449 | 1.4220469 | 0.5576784 |
| transit_intensity | 3.5481285 | 8.6480176 | 0.4102823 |
| co2_capita | 0.8639931 | 0.6056498 | 1.4265556 |
palma <- cbind(rownames(palma), palma)
write.xlsx(palma, "palma.xlsx", overwrite = TRUE)
promedios <- data.frame(dist_capita = sum(UTAM$dist)/sum(UTAM$pop),
walk_intensity = sum(UTAM$dist_walk)/sum(UTAM$pop),
auto_intensity = sum(UTAM$dist_auto)/sum(UTAM$pop),
transit_intensity = sum(UTAM$dist_transit)/sum(UTAM$pop),
ges_capita = 1000*sum(UTAM$co2)/sum(UTAM$pop))
kable(promedios, caption = "Valores promediados de los indicadores", format.args = list(big.mark = ","))
| dist_capita | walk_intensity | auto_intensity | transit_intensity | ges_capita |
|---|---|---|---|---|
| 12.71442 | 1.161807 | 1.72813 | 7.015431 | 0.6600684 |